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Abstract 

In this contribution, the unsteady, viscous and two 
dimensional flow over a plunging NACA0012 airfoil was 
simulated using a high-order numerical method associated 
with a moving frame of reference technique. The goal of 
using a high-order numerical method is to characterize 
heaving airfoil wakes with very fine temporal and spatial 
resolutions. This characterization is essential in 
understanding the transition process experienced by 
oscillating airfoils and which can enable more control on the 
transition. The response of motionless and periodically 
forced airfoils was considered separately. The flow was 
investigated at a mean incidence a = 0°, Reynolds number Re 
= 10 4 , and over a range of heaving frequencies and 
amplitudes. The resulting airfoil wakes were analyzed with 
respect to their frequency response and type. It is shown that 
for forced airfoils three regimes of frequency response exist: 
a natural regime, a lock-in regime, and a harmonic regime. 
The transition between these regimes was presented at 
constant frequency and constant amplitude of the forcing 
oscillations. The frequency regimes were then explored in 
detail and related to the wake types and to the wakes 
transition. The results obtained were in accordance with 
previously published experimental and computational 
investigations. Furthermore, a more detailed wake 
classification was presented, where two more wake 
configurations were captured. The role of the trailing- and 
leading-edges in the creation of multiple-vortices-per-half- 
cycle shedding mode was investigated. Moreover, the 
interaction of the natural shedding frequency of the airfoil 
with the imposed forcing frequency and their contribution in 
the creation of multiple-vortices-per-half-cycle shedding 
mode was also investigated. These results showed that 
frequency interactions can be related to the appearance of 
the multiple-vortices-per-half-cycle modes. 
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Introduction 

Oscillating airfoils have received much attention and 


interest in different contexts and applications due to 
the unsteadiness and the complex nature of these 
flows, rendering their numerical simulation very 
challenging. Unsteadiness and its consequences are 
fundamentally important problems which concern a 
broad range of applications. These undesirable 
consequences can range from flutter (Scalan 1971), to 
vortex-induced-vibrations (VIV) (Williamson 2004), 
and flow-induced- vibrations (FIV) (Belvin 1977; 
Naudascher 1994). 

Oscillating airfoils include: heaving (or plunging), 
pitching and flapping airfoils. Purely heaving airfoils 
have received less attention compared to pitching 
airfoils, which is mainly due to the applications 
involving dynamic stall on helicopters which are more 
concerned by the pitching motion. Nevertheless, 
heaving airfoils exhibit a rich variety of wakes which 
are classified into three types as drag-producing, 
neutral and thrust-producing wakes (Freymuth 1988; 
Lai 1999; Young 2005; Medjroubi 2011b; Medjroubi 
2012), depending on the frequency and amplitude of 
the heaving oscillations. At low reduced frequencies, 
the wake is a Karman vortex street, and the wake 
produces drag; while at high reduced frequencies, the 
wake becomes a reversed Karman vortex street. 
During the transition from drag to thrust-producing 
wakes, a neutral mode appears where the produced 
thrust equals the inherent drag. In some cases, this 
mode is characterized by the shedding of more than 
one vortex pair per oscillation cycle. This shedding is 
referred to as the multiple-vortex-per-half-cycle mode. 
The details of the wake structure formation, and the 
lift and drag production are not well understood, 
especially the multiple-vortex-per-half-cycle mode 
(Young 2001; Young 2005). The multiple-vortex-per- 
half-cycle mode origin is thought to be the result of the 
interaction between the natural shedding of the airfoil 
and the imposed heaving frequency (Young 2005; 
Medjroubi 2012). 

The forced wake of the airfoil exhibits also three 
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different frequency regimes. These regimes depend on 
the forcing frequency (and the forcing amplitude) and 
the relation between the forcing frequency and the 
natural shedding frequency of the airfoil. The natural 
frequency (or Strouhal frequency) is defined as the 
frequency of shedding of vortices of the unforced 
airfoil (Young 2005; Karniadakis 1989). The three 
regimes exhibited are called natural, harmonic, and 
lock-in regimes, and have been extensively studied 
and observed for oscillating cylinders (Ongoren 1988a; 
Ongoren 1988b; Karniadakisl989; Cheng 2001; Young 
2004), but less explored for oscillating airfoils (Young 
2005; Young 2004). The objective of this paper is to 
explore these regimes and to determine the relation 
between the frequency regimes and the wake types. 

Spectral/hp Element Methods was used in this work 
combined with a moving-frame of reference technique 
to investigate the problem of the flow around an 
oscillating airfoil (Medjroubi 2011b). All the 
simulations presented in this work are Direct 
Numerical Simulations (DNS), in order to capture all 
the flow details and particularly the spatio-temporal 
flow evolution. To the authors knowledge, it is the 
first time that the frequency-regimes were investigated 
for a heaving airfoil using a high-order numerical 
approach, as most of the studies for this type of 
motion were concerned with the wake configuration 
(Young 2005; Anderson 1998; Ellenrieder 2003). 

In this paper, the incompressible, 2D and unsteady 
flow over a stationary and a heaving NACA0012 
airfoil was simulated. The Reynolds numbers 
considered are moderate (Re=10 4 -3xl0 4 ), with Re 

number based on the airfoil chord length). A range of 
heaving amplitudes and frequencies was simulated in 
order to expose the frequency regimes introduced 
earlier and to achieve a transition between these 
regimes. The dependency of the frequency regimes on 
the relation between the natural and forcing frequency 
was explored as the origin of the multiple-vortex-per- 
half-cycle shedding mode. The Spectral /hp Element 
Method and the moving-frame of reference technique 
are introduced in Section 2. The characterization of the 
wake modes in function of the interaction of the 
natural shedding frequency and the imposed heaving 
frequency interaction are presented in Section 3. 
Finally, Section 4 containes a summary of the findings 
of this paper. 

The Numerical Method 

The main idea behind using high-order numerical 
methods is to achieve high accuracy, resolving all the 
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flow length scales at a reduced cost and avoiding the 
extra costs of remeshing or excessive grid resolution 
(Nastase 2006). Unlike pure spectral methods and 
finite elements methods, Spectral/hp ElementMethods 
(SEM) allows the use of two refinement techniques 
simultaneously: the h-refinement, denoting the 

increase of the number of elements, and the p- 
refinement, denoting the increase of the polynomial 
order of the approximation (Karniadakis 2005). The p- 
refinement is a powerful tool, as it offers the 
possibility of increasing the simulations accuracy 
without any remeshing by increasing the order of the 
approximation in all or parts of the computational 
domain, which can be considered as a way of 
implementing adaptivity. For more details about the 
derivation and applications of SEM, the reader can 
refer to (Karniadakis 2005; Karniadakis 1991; 
Medjroubi 2011a). To implement the airfoil motion, the 
Navier-Stokes equations and the boundary conditions 
were written in a moving frame of reference. In the 
following, the principles of the Spectral/hp Element 
Method and the moving frame of reference are 
presented. Readers interested in the details of the 
derivation and the implementation of the method can 
refer to the monograph by S. Sherwin and G. Em. 
Karniadakis (Karniadakis 2005). 

Spatial and Spectral Discretization 

The 2D Navier-Stokes equations for an incompressible, 
unsteady and viscous flow are written as: 

V u = 0, (1) 

— + (u-V)u =-^- + v V 2 u + F. (2) 

dt y ’ p 

where: u is the velocity vector, p the pressure, v the 
kinematic viscosity and F are the body forces. 

The physical domain is divided into triangular and/or 
quadrilateral (tetrahedral and/or hexahedrals) sub- 
domains for 2D (3D) configurations. A spectral 
expansion in the form of Jacobi polynomials of mixed 
weight and of order p, was used to represent the 
solution variable within each sub-domain. 

Polynomial Bases/Expansions 

For a polynomial expansion, a basis u mn (r,s) needs to 
be defined in order to approximate the function / (x, y) 
over i triangular sub-domains using a C° continuous 
expansion of the form (for basis used for rectangular 
sub-domains refer to (Karniadakis 2005)): 

f(x,y) = T'L'L U m,n U m,n( r ( X ’y'>’ S ( X ’yV ( 3 ) 

i m n 


86 



Review of Applied Physics Vol. 2 Iss. 4, December 2013 

where: u l mn are the expansion coefficients in the i th 
sub-domain, (x,y) are the spatial coordinates and (r,s) 
are the local coordinates within the sub-domains. N is 
the total number of sub-domains, N l9 N 2 are the 
number of the quadrature points used in the r ands 
directions; respectively. The space spanned by the 
local coordinates system is defined as follows: 

]} ={( r>s ) | -l<r,s ; r + s<0} (4) 

where the orthogonal expansion bases used are 
Dubiner's modified bases, defined as: 

( 5 ) 

pa^(x) is the n th -order Jacobi polynomial in the [-1,1] 
interval, which satisfies the orthogonality relationship: 

} PZ-fi (x)P^ (x)(l - x) a (1 + xf dx = S : , (6) 

-1 

where: 8™ is the Kronecker-delta. 

The bases are decomposed into boundary and interior 
modes, thus allowing the construction of a global C° 
expansion. Integration and differentiation were 
performed at element level, and an element-mapping 
which allows the generalization of the local operations 
in a standard region to elements of general shapes is 
defined. To extend these techniques to a C° multi- 
dimensional basis, global operations as matrix 
numbering, connectivity and assembly were 
introduced. This spatial discretization thereby 
introduced is independent from the Navier-Stokes 
solver. The solver used in this work will be introduced 
in the next section. 

Navier-Stokes Solver: High Order Splitting Scheme 

The temporal discretization of the Navier-Stokes 
equations was achieved via a splitting scheme. This 
temporal scheme requires three steps to determine the 
fields for the next time step using the velocity at the 
previous time step. The linear terms in the Navier- 
Stokes equation are treated implicitly, while the 
nonlinear terms are treated explicitly. A Poisson 
equation is solved for the pressure, and this is 
achieved by forcing the incompressibility constraint 
over an intermediate velocity field. The final velocity 
field is obtainedby solving a Helmholtz equation. For 
more details about this scheme, refer to (Karniadakis 
1991). 

The Dimensionless Navier-Stokes Equations 

The non-dimensional Navier-Stokes equations were 
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obtained by making use of the 
dimensionless variables: 

following 

* x * y 

X =c' ? =? 

( 7 ) 

t* _ £H® u * _ JL 

c U 0 

( 8 ) 

* V * p 

V = rr ’ P = 2 ’ 

( 9 ) 

U 0 pu\ 

where: U 0 is the free-stream velocity, c is 

the chord 

length of the airfoil, u and v the horizontal and 
vertical velocity components. The Reynolds number is 

defined as Re = . 

V 


After introducing the non-dimensional variables, the 
dimensionless Navier-Stokes equations are obtained: 

o' 

II 

* 

3 

> 

(10) 

dll / * 1 * 

bill -Vlu =-V» + — V 2 u . 

dt ' > Re 

(11) 


For simplification the * sign is dropped. 
Airfoil Motion 


To implement the oscillatory motion of the airfoil, the 
Navier-Stokes equations are written in a moving frame 
of reference and an equation is derived for the airfoil 
motion. The resulting system is then solved. The 
method is shortly described here, and for further 
details, refer to (Li 2002). 

Consider that the airfoil performs a translational 
motion d{t) defined in the fixed or absolute frame of 
reference (X a , Y a ) : 

d = (a(x),b(y)). ( 12 ) 

A moving frame of reference is then attached to the 
airfoil, (x m ,y m ) which can be definedas a function of 
the absolute frame as: 

X a =a(x) + x m , (13) 

Y a - b(x) + y m . (14) 



FIG. 1 SCHEMATIC REPRESENTATION OF THE AIRFOIL FOR 
THE HEAVE MOTION, WHERE THE VERTICAL 
DISPLACEMENT y(t) VARIES SINUSOIDALLY. 0 IS THE MEAN 
INCIDENCE (INITIAL ANGLE OF ATACK), U x IS THE FREE- 
STREAM VELOCITY AND c THE CHORD LENGTH. 


87 


ww w. seipub . or g / rap 

Using the definition of the coordinates (x m ,y m ) given 
by the equations (13) and (14), the Navier-Stokes 
equations are then re-written in the moving frame of 
reference as: 



< 

c 

II 

© 

(15) 

du , 

— + (u 
dt v 

v 1 O 

•V)u = -V/? + — v u + Q(u,t), 
Re 

(16) 


Q(u,t) = -A t d, 

(17) 


where: A T d results from the unsteady translational 
motion. The boundary conditions are also written in 
the moving frame of reference, and expressions are 
derived for Neumann and Dirichlet boundary 
conditions. For the present simulations, the airfoil is 
considered as rigid and forced to oscillate in heave 
motion as illustrated in Figure 1. The heave motion is 
defined as: 

d (i t ) - y(t ) = y 0 + a cos(2 n ft). (18) 

where: y(t) is the time dependent vertical motion, y 0 
is the initial heaving amplitude, / = f* = fc / and a 
the dimensionless frequency and amplitude of 
oscillations; respectively. 

The airfoil is set at a constant initial angle of attack a 
and forced to oscillate vertically in a sinusoidal fashion. 
The boundary conditions at the domain surfaces are 
inflow conditions at the entry of the domain and 
outflow conditions at the domain exit. Periodic 
boundary conditions were used for the upper and 
lower part of the domain. The boundary conditions at 
the airfoil surface are no-slip conditions. 

Flow Solver 

The general purpose Navier-Stokes solver NeKxar was 
used, based on the Spectral /hp Element Method. The 
solver was validated for many flow geometries, 
among those cylindrical (Carmo 2008) and rectangular 
(Li 2002). The solver was validated for the airfoil 
geometry for steady and unsteady flows (Medjroubi 
2001b; Medjroubi 2011a; Medjroubi 2012) and the 
results obtained were in agreement with previously 
published numerical results (Shatalov 2006; Lian 2007; 
Gopalan 2008) was found. All the simulations of 
oscillating airfoils were started from fully converged 
motionless airfoil solutions. 

Results 

The results presented here were obtained using a 2D 
unstructured grid composed of 4220 elements. Details 
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about the numerical grid, the boundary conditions 
used and the parameters of the SEM method can be 
found in (Medjroubi 2011b). 

Unforced Airfoil Wakes 

The flow is simulated at a = 0° and at a Reynolds 
numbers Re = 10 4 . The flow exhibits a typical Karman 
vortex-street configuration, attached over all the 
surface of the airfoil and leaves the airfoil parallel to 
the tr ailing-edge. The vortices are shed at a certain 
distance from the trailing-edge due to the wake 
instability (Schnipper 2009) and this distance decreases 
as the Reynolds number is increased. The wake is a 
drag-producing wake, and this is confirmed by the 
value of the mean thrust coefficient. 

A Fast Fourier Transform (FFT) was performed on the 
time series u(t ) at Re = 10 4 and a = 0° in order to obtain 
its frequency content. This technique was used to 
characterize the frequency content of the airfoil 
response either for the motionless airfoil or for the 
heaving airfoil. Fig. 2 illustrates the FFT of the time 
trace of the horizontal velocity u(t ) . The peak 
frequency k - k nat - 7.85 is called the natural shedding 
frequency, or frequency at which the vortices are shed 
into the wake of the airfoil. It characterizes the 
response of the airfoil to the incoming inflow at a 
certain Reynolds number (Green 1995). The frequency 
response contains also some other weaker peaks at the 
superharmonics of the natural shedding frequency. 
This behavior is typical of a motionless object response 
to an incoming inflow (Blevins 1977; Karniadakis 1989; 
Williamson 1995; Mittal 1999). 

Forced Airfoil Wakes 

The response of the airfoil wake to a harmonic forcing 
was presented in this Section. The airfoil was forced to 
oscillate in heave at a prescribed frequency and 
amplitude at Re = 10 4 and a mean angle of attack a = 
0°. The imposed heaving frequency was labeled 
forcing frequency,^ . The response frequency of the 

flow to forcing was obtained by analyzing the 
frequency content of the temporal signal of the 
horizontal velocity u(t ) using a FFT. The dominant 
frequency in the response signal was labeled k r . The 
frequency and amplitude values considered in this 
simulation were chosen to exhibit the heaving airfoil 
wake in different frequency-response regimes. 

These regimes are the natural regime, where the 
response frequency is equal to the natural shedding 
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frequency of the motionless airfoil. The lock-in regime 
corresponds to the response frequency locked to the 
forcing frequency. Finally, the harmonic regime 
corresponds to the response frequency being a mixture 
of the natural and the forcing frequency. 

Constant Frequency Regime Transition 

The forcing frequency was kept constant and the 
forcing amplitude was varied. The frequency content 
of the horizontal velocity time signal was extracted at 
a location in the near wake (x= 1.63, y= 0.09) and its 
frequency content was then determined by performing 
an FFT. This location was chosen according to the 
assumption that it is in the near wake that the forcing 
and natural frequency compete with each other. The 
result of this frequency competition or interaction 
determined the final state of the wake (Karniadakis 
1989). 



FIG. 2 FAST FOURIER TRANSFORM OF u(t) FOR A 
MOTIONLESS AIRFOIL AT Re = 10 4 AND a = 0°. 

The time trace of the horizontal velocity and the FFT 
analysis of the natural regime (or mode) are shown in 
Fig. 3a and Fig. 3b. This regime was obtained at 
k — 1.96 (corresponding to k / k nat = 1 / 4 ) and h =0.001. 
The dominant peak occurred at a frequency equal to 
the natural shedding frequency k nat =7.85 and its first 
super- harmonic. Another peak occurred at the forcing 
frequency but weaker than the ones at k nat and its 
super-harmonic. Some other weak peaks were present 
around the dominant peaks, and this behavior was 
also observed for a cylinder in a non-lock-in state and 
reported in (Karniadakis 1989). Increasing the forcing 
amplitude to h= 0.04 and keeping the frequency of 
oscillation constant have resulted in a dramatic change 
in the frequency response of the flow. Fig. 3c shows 
the time-series of u and Fig. 3d shows its FFT analysis. 

The peak frequency is k r = 5.93 corresponding to a 
linear combination of the natural shedding frequency 
and the forcing frequency ( k r = k nat /2 + kj- ), which 


indicates a harmonic regime where k nat and kj 
inter acte together. The result of this interaction is that 
the peak frequency is some mixture of the natural 
shedding frequency and the imposed forcing 
frequency. Other weaker peaks are present at some 
other combinations of these two frequencies ( k nat + 
1/4 k nat + ky). 


Compared to the frequency response of the natural 
regime, it is important to note that the k r peak 
becomes more important and the peak corresponding 
to the natural shedding frequency k nat is less 
prominent, indicating that the forcing frequency gains 
importance at the expense of the natural shedding 
frequency. 
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FIG. 3(A) TIME TRACE OF THE HORIZONTAL VELOCITY u 
AND ITS (B) FAST FOURIER TRANSFORM h =0.001, k = 1.96. (C) 
TIME TRACE OF u AND ITS (D) FAST FOURIER TRANSFORM 
AT h = 0.04, k = 1.96. (E) TIME TRACE OF u AND ITS (F) FAST 
FOURIER TRANSFORM AT h =0.1, k = 1.96. ALL SIMULATIONS 
ARE AT RE = 10 4 , a = 0°, AND u(t) IS PLOTTED AT THE 
LOCATION (x = 1.63, y = 0.09). 

A further increase in the forcing amplitude to h =0.1 
results in a lock-in regime where the response 
frequency is locked to the forcing frequency. Fig. 3e 
illustrates the time trace of the horizontal velocity u . 
Fig. 3f shows the FFT of the time trace of u where the 
peak frequency now corresponds to the forcing 
frequency kf = k r = 1.96 .There are less important 
peaks at the super-harmonics of the forcing frequency. 
This is typical for a lock-in regime, as this regime is 
similar to the non-forced case in its frequency response 
(see Fig. 2) (Karniadakis 1989). 

The frequency regime transition described above 
embodies another transition in the wake type (or 
pattern). This transition is from a drag-producing 
wake to a thrust- producing wake via a neutral wake 
at constant frequency. Fig. 4 shows the vorticity 
distributions of the simulations described above. 
These wake types are typical of heaving airfoils and 
the transition between them has been reported in the 
literature (Lai 1999; Young 2004; Young 2005; Platzer 
2008; Medjroubi 2011b; Medjroubi 2012). In Fig. 4a, the 
wake exhibits a Karman vortex-street, where a pair of 
counter-rotating vortices is shed peroscillation cycle. 
This is a drag-producing wake. The wake type 
changes as the amplitude of the heaving oscillation 
increases (see Fig. 4c). This wake type is called a 
neutral wake, as the inherent drag balances the 


produced thrust. 

This wake is also labeled multiple-vortices-per-half- 
cycle mode, as more than one vortex-pair is shed per 
oscillation cycle. This mode is not well understood, 
and the mechanism by which more than a vortex is 
shed per half-cycle is not well known (Young 2005). 


It was proposed here that this mode (or wake type) 
results from the interaction between the shedding 
frequency and the imposed frequency, as already 
explained earlier in this Section by means of FFT 
analysis. This result will be further confirmed by the 
simulations of the frequency transition at constant 
amplitude in the next Section. This explanation for the 
multiple-vortices-per-half-cycle regime was also 
provided by Young (Young 2005) who briefly 
discussed this mode and did not provide any 
frequency analysis, and then concluded only based on 
observation of the time trace of u for one cycle of 
oscillations. A further increase of the forcing 
amplitude results in a transition to a thrust- producing 
wake (see Fig. 4e). Note here that more than one pair 
of counter-rotating vortices is shed per oscillation- 
cycle. This thrust-producing wake is different from the 
typical thrust-producing wake, as the latter looks like 
the drag-producing wake in Fig. 4a. This new thrust- 
producing wake is due to the increase in the 
amplitude of the heaving oscillation. The increase in 
the amplitude triggers the shedding of more vortices 
per oscillation cycle (Young 2004). This is mainly due 
to the fact that at higher amplitudes the leading-edge 
also contributes to the vortices shed at the trailing- 
edge. This is confirmed by the presence of vortices on 
the upper and lower airfoil surface. Note that for the 
multiple-vortices-per-half-cycle for the neutral regime 
(see Fig. 4c) there are no contributions from the 
leading-edge, although the flow is slightly detached 
close to the trailing-edge region. 
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(f) 


FIG. 4 VORTICITY DISTRIBUTION AT A = 0°, RE = 104 FOR (A)A 
DRAG-PRODUCING WAKE ( k = 1.96, h = 0.001),(B) VIEW OF THE 
NEAR WAKE REGION.(C) A NEUTRAL WAKE ( k = 1.96, h = 
0.04),(D) VIEW OF THE NEAR WAKE REGION. (E) A THRUST- 
PRODUCING WAKE ( k = 1.96, h = 0.1), (F) VIEW OF THE NEAR 
WAKE REGION. 
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k=3.92 
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(c) 



Reduced frequency k 

(d) 

FIG. 5 (a) TIME TRACE OF THE HORIZONTAL VELOCITY u 
AND ITS (b) FAST FOURIER TRANSFORM AT k = 1.96, h = 0.05. (c) 
TIME TRACE OF THE HORIZONTAL VELOCITY u AND ITS (d) 
FAST FOURIER TRANSFORM AT k = 3.92, h = 0.05. ALL 
SIMULATIONS ARE AT RE =10 4 , a = 0°, AND u{t) IS PLOTTED AT 
THE LOCATION ( X = 1.63, y = 0.09). 


Constant Amplitude Regime Transition 

In this Section, the forcing amplitude is kept constant 
at h = 0.05 and the forcing frequency is increased. The 
frequency content of u(t ) is analyzed at the same 
location as for the constant frequency regime 
transition. 

The time trace of u at h = 0.05 and k = 1.96 is 
illustrated in Fig. 5a. The peak frequency of the time 
trace of u corresponds to a mixture of the natural 
shedding frequencyand the forcing frequency 
( k r = 1/4 k nat + kj = 3.92 ). A second frequency peak 

corresponds to another combination of these 
frequencies ( k r = 1/2 k nat + kf = 5.88 )(see Fig. 5b). 

There are weaker peaks at the natural shedding 
frequency and the forcing frequency and their sum 
( k = k nat +kf = 9.81 ). This response is typical of a 

harmonic regime, as shown earlier in this paper for the 
simulation at ( k = 1.96, h = 0.04). However, it is 
important to note that for the present simulation ( k = 
1.96, h = 0.05) the peak at the natural frequency 
k nat =7.85 is more important than the one at the 
forcing frequency^ =1.96, while this is inverted for 

the simulation at ( k = 1.96, h = 0.05). This may indicate 
that the response of the airfoil to the forcing motion 
depends on the amplitude of the oscillations. 

Increasing the forcing frequency to k y =3.92 results in 

a transition to a lock- in regime (see Fig. 5d), where the 
dominant frequency peak is equal to the forcing 
frequency^ =3.92 . The time trace of the velocity is 
presented in Fig. 5c. During the transition from the 
harmonic to the lock-in regime, a wake transition 
occurs from a neutral to a thrust producing wake. The 
vorticity contours of the neutral wake (not presented 
here) show that the wake presents a multiple-vortices- 
per-half-cycle regime. As for the constant frequency 
transition, the neutral multiple-vortices-per-half- cycle 
regime has no vortices contribution from the leading- 
edge. The flow is attached over all the airfoil surface. 
The apparition of this mode in the harmonic region, in 
both simulations at constant forcing frequency and 
constant forcing amplitude, confirms that the origin of 
this mode is related to the natural and forcing 
frequencies interaction. Increasing the frequency of 
oscillations results in a wake transition to a thrust- 
producing wake (not presented here). In the wake, one 
pair of counter-rotating vortices is shed per oscillation 
cycle. The wake configuration is typical of a thrust- 
producing wake (Young 2001; Young 2005), and 
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although there is a slight contribution of leading-edge 
to the shed vortices, this does not result in a multiple- 
vortices- per- half- cycle configuration. This 
phenomenon is discussed in the next Section. 

Frequency Regimes Classification - a Frequency 
Pattern? 

Several simulations were conducted at different 
combinations of frequency and amplitude of the 
forcing oscillations in order to obtain a frequency 
regime classification. The simulations were conducted 
at low heaving amplitudes and at moderate to high 
heaving frequencies. This region of the k / k nat vs h plot 
is associated with aero-elastic and aeroacoustic 
phenomena and has been not well explored in the 
literature. The classification was obtained by 
identifying the different frequency regimes using the 
FFT analysis and expressed in function of the 
amplitude and the oscillation frequency normalized 
by the natural frequency of vortex shedding, k / k nat 
(see Fig. 6). 

At low forcing frequencies ( h < 10 -2 ), an increase in the 
frequency at a constant amplitude triggers a transition 
from the natural regime to the lock-in regime without 
going through a harmonic regime. This is also valid 
for a transition varying the amplitude at a constant 
forcing frequency for k / k nat > 0.5. This is because 
above k / k nat = 0.5, any harmonic of the forcing 
frequency will be higher than the natural shedding 
frequency and the flow will lock to the forcing 
frequency. This is mainly due to the trailing-edge of 
the airfoil which dictates the separation point for 
k / k nat >0.5 and enforces the dominance of the forcing 
frequency (lock-in) even in the presence of leading- 
edge vortex shedding, as it is observed in the case of 
the constant amplitude regime transition described in 
Section 3.2.2. The vorticity contours at h = 0.05 and k = 
3.92 (corresponding to k / k nat = 0.5, represented on Fig. 
6 with a star pointed to by a red arrow), not presented 
here, show that vortices are shed at the trailing-edge 
when the airfoil is at its maximum y max and minimum 
y min vertical position. This is an indication that the 
trailing-edge motion imposes the shedding process, 
even if vortices sare present at both the upper and 
lower airfoil surfaces. 

Fig. 7 shows the vorticity contours of the wake at h = 
0.1 and k = 1.96 (k/k nat = 0.25, represented by a black 
circle pointed to by a red arrow on Fig. 6). The 
shedding of vortices occurs not only when the airfoil is 
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atjmax or Tmin but also when the airfoil is at y = 0. This 
explains how more than one vortex-pair is shed per 
oscillation cycle. In this regime the trailing-edge 
heaving motion does not control the shedding process. 
It was believed that the interaction between the 
natural shedding frequency and the imposed 
frequency is responsible for this particular shedding 
regime. The response frequency in this regime that is 
k = 1.24, can be written as a combination of both k na1 
and kf . On the other hand, and when the streamlines 

were plotted around the trailing-edge region (see Fig. 
8), a recirculation region was observed on both sides 
of the trailing-edge during the whole heaving cycle. 
This circulation region varies in strength and position 
during the heaving cycle, thus creating a blunt-like 
body configuration at the trailing-edge. The variation 
of this bluntness with the heaving cycle has a 
frequency associated with it. This frequency be neither 
the natural nor the forcing frequency but a mixture of 
both. 



FIG. 6 CLASSIFICATION OF THE WAKE-TYPES AS A FUNCTION 
OF THE AMPLITUDE h AND THE NORMALIZED FORCING 
FREQUENCY k / k nat . THE TRIANGLES ARE FOR DRAG- 
PRODUCING WAKES, WITH ONE VORTEX PAIR SHED PER 
OSCILLATION CYCLE. THE STARS REPRESENT THRUST- 
PRODUCING WAKES WITH ONE VORTEX PAIR SHED PER 
OSCILLATION CYCLE. THE SQUARES ARE FOR NEUTRAL 
WAKES WITH MULTIPLE-VORTICES-PER-CYCLE SHEDDING. 
THE OPEN CIRCLES REPRESENT NEUTRAL WAKES WITH ONE 
VORTEX PAIR SHED PER OSCILLATION CYCLE. FINALLY, THE 
BLACK CIRCLES ARE FOR THRUST-PRODUCING WAKES WITH 
MULTIPLE-VORTICES-PER-CYCLE SHEDDING. THE 
FREQUENCY REGIMES BOUNDARIES ARE APPROXIMATELY 
INDICATED BASED ON THE SIMULATIONS PERFORMED IN 
THIS PAPER. 
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(c) (d) 


FIG. 7 VORTICITY CONTOURS AT a = 0°, Re = 10 4 . THE 
HEAVING AMPLITUDE IS h = 0.1 AND THE FREQUENCY IS k = 
1.96. THE FIGURES FROM (A) TO (D) REPRESENT ONE 
HEAVING CYCLE, AND Y IS THE VERTICAL POSITION OF THE 
TRAILING-EDGE. (A) Y=0, (B) Y = j; min , (C) Y=0, (D) Y = j; max . 

In the region defined by h > 10 -2 and k / k nat < 0.5, the 
transition at constant frequency or amplitude from the 
natural to the lock-in regime passes via the harmonic 
regime, as observed in the previous Section. 

The classification of the different wakes in terms of the 
frequency regimes provides an idea of the complexity 
of the frequency regimes transition and the 
frequencies interplay in the case of a heaving airfoil. 
There is a rich variety of wake-types and frequency 
regimes. It can be concluded that the wake-type is 
strongly related to the frequency regime. For example, 
wake-types with multiple-vortices-per-half-cycle of 
oscillation are in a region where the force and the 
natural frequency of vortex shedding interact together 
(the harmonic regime). Wake-types with one vortex- 
pair shed per cycle are in the region where one distinct 
frequency is in control. 

Note also that the phase diagram is not symmetric, as 
it is the case for the cylinder (Karniadakis 1989; Young 
2005). This can be attributed to the presence of the 
sharp trailing-edge of the airfoil. 
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FIG. 8 STREAMLINES AT a = 0°, Re = 10 4 , h = 0.1 AND k = 1.96 
DURING ONE OSCILLATION CYCLE. (A) Y=0, (B) Y = j; min , (C) 
Y=0, (D)Y = j, max ,(E)Y=0. 

The wake classification presented in this paper agrees 
very well with the one presented by J. Young in 
(Young 2005) at Re = 2 x 10 4 , although the author 
noted the presence of only three wake types in his 
diagram and no frequency analysis was presented to 
confirm the nature of the frequency regimes. A less 
similar diagram for the pitching airfoil was shown by 
T. Schnipper et al. (see (Schnipper 2009), Fig.2). 
Although it is important to note that for both the 
pitching and heaving cases, there exist a region where 
the shedding of more than a vortex pair per oscillation 
occurs and that this region is confined to high 
amplitudes and low oscillation frequencies. 

Conclusion 

The simulations of two dimensional unsteady flow on 
both a motionless and a heaving airfoil using a high- 
order spectral CFD method were successful. The high 
spatial and temporal resolution achieved by the 
Spectral /hp Element Method permitted a detailed 
analysis of the flowover the airfoil and the near wake. 
The high-order numerical methods enabled the 
simulation and the characterization of the frequency 
response regimes. These regimes are: (a) the natural 
regime, (b) the lock-in regime, and (c) the harmonic 
regime. These frequency-response regimes are related 
to the shedding process through which the wake 
undergoes a transition from a Karman street (a drag- 
producing wake) to a reversed Karman street (a 
thrust-producing wake). The transition between the 
different frequency regimes was simulated at both 
constant frequency and constant amplitude. This 
transition is successfully related to the wake-types 
regimes and the transition between them. It has been 
found that the frequency regimes are strongly related 
to the wake type exhibited. Wake-types with multiple- 
vortices-per-half-cycle of oscillation were found in 
harmonic regimes and wake-types with one vortex- 
pair shed per cycle were in the region where one 
distinct frequency was in control (lock-in and natural 
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regimes). Regimes where the leading-edge vortices 
contribute to the shedding process were also 
simulated. 
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